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Abstract 


1 Motivation: 

This paper presents libRoadRunner, an extensible, high-performance, cross-platform, open-source soft¬ 
ware library for the simulation and analysis of models expressed using Systems Biology Markup Language 
(SBML). SBML is the most widely used standard for representing dynamic networks, especially biochem¬ 
ical networks. libRoadRunner supports solution of both large models and multiple replicas of a single 
model on desktop, mobile and cluster computers. 


2 Results: 

libRoadRunner is a self-contained library, able to run both as a component inside other tools via its C++ 
and C bindings andnteractively through its Python interface. The Python Application Programming 
Interface (API) is similar to the APIs of Matlab and SciPy, making it fast and easy to learn, even for 
new users. libRoadRunner uses a custom Just-In-Time (JIT) compiler built on the widely-used LLVM 
JIT compiler framework to compile SBML-specified models directly into very fast native machine code 
for a variety of processors, making it appropriate for solving very large models or multiple replicas of 
smaller models. libRoadRunner is flexible, supporting the bulk of the SBML specification (except for 
delay and nonlinear algebraic equations) and several of its extensions. It offers multiple deterministic 
and stochastic integrators, as well as tools for steady-state, stability analyses and flux balance analysis. 


3 Availability and Implementation: 

We regularly update libRoadRunner binary distributions for Mac OS X, Linux and Windows and license 
them under Apache License Version 2.0. http://www.libroadrunner.org provides online documentation, 
full build instructions, binaries and a git source repository. 
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5 Introduction 


Dynamic network models (Saurol , 2014h of metabolic, gene regulatory, protein signaling and electropliysio- 
logical models require the specification of compon ents, in teract io ns, compartments and kinetic parameters. 
The Systems Biology Markup Language (SBML) ( Hucka et all 12004 ) has become the de facto standard for 
declarative specification of these types of model (see SBML.org and refs.). 

Po pular tools for the development, simulation and analysis of models specified in SBM L include CO- 
PASI ( Hoops et al. I . 2006h . Systems Biology W orkbench (SBW) (Bergmann and SauroL 2006), T he Systems 
Biol ogy Simulation Core A lgorithm ( TSBSC) ( Keller et all [20 1 .'ll ) . Jarnac (Sauro and Fell . 2000t) . libSBML- 
Sim (Takizawa et a,L . 201B), SOSLib (SOSLib Team . 20l4 ). iBioSim ( Mvers et oil 2009 ). PySCeS ( Olivier et all 


2005), and VirtualCell ( Moraru et all 20081) . Some of these applications are stand-alone packages designed 
for interactive use, with limited reusability as components in other applications. Very few are flexible 
component libraries. Currently, none are fast enough to support emerging applications t hat req uire l arge - 


scale simulation of network dynamics. For example, multicell virtual-tissue simulations ( Hester et al.. 120111) 


often require simultaneous simulation of tens of thousands of replicas of models residing in their cell ob¬ 
jects and interacting between cells. In addition, optimization methods require generation of time-series 
for tens o f thousands of replicas to explore the high-dimensional parameter spaces typical of biochemical 


networks (Bout eille r et al, , 2014). 


We designed libRoadRunner to provide: 1) Efficient time-series generation and analysis of large or mul¬ 
tiple SBML-based models; 2) A comprehensive and logical API; 3) Interactive simulations in the style of 
IPython and MATLAB; and 4) Extensiblity. 

Most existing SBML simulation engines use built-in interpreters to parse and execute SMBL model spec¬ 
ifications. Interpreted execution is simple and flexible, but much slower than execution of compiled code. 
Other simulation engines generate compiled executables from SBML by first converting SBML-specified mod¬ 
els into a general-purpose-language representation. The engines then call an external compiler to translate 
the general-pur pose-language into an execu table shared library to load at run tim e. E.q., roadRunner in 
the SBW suite ( Bergmann and SauroL 20061) converts SBML into C# (see § 1.4 of ( Alfred V. Ahol . 19861) ), 
then compiles the Cff using a compiler from the .NET distribution. This approach generates relatively fast 
executables. However, it requires distribution of a separate compiler or a redistributable runtime, reducing 
portability. 

A more efficient approach to SBML-to-executable compilation uses a specialized just in time JIT com¬ 
piler, to compile SBML into an optimized Intermediate Language ( IL ) r epresentatio n and the IL code into 


native executable machine code directly in-memory. Ackermann et al. ( Ackermann et al. . 1 2009 ) used JIT 


compilation to generate CUDA code from SBML and execute it on an Nvidia GPU. libRoadRunner sup¬ 
ports execution of a broad range of_SBML models on CPUs using a custom-built JIT compiler (based on the 
LLVM JIT compiler framework ( Lattner and Advil 2004 )) which translates SBML into highly-optimized 
executable code for a broad range of processors. LLVM-based compilers are small, so all JIT operations 
occur in memory, without external file or compiler access, ensuring fast, self-contained simulations and a 
relatively small distribution package. 


Capabilities libRoadRunner supports time-course simulation or deterministic and stochastic models. In 
addition it supports steady state, stability analysis and flux balance analysis of SB ML-specified models. 
For flux balance analysis libRoadRunner support the FBC SBML extension (jOlivier and Bergmanm 2012). 
It supports almost the entire SBML L3V1 specification, including composite model composition and the 
distribution package. Its only non-compliance is its lack of support for delay equations and non-linear 
algebraic rules. 
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Portability Because new hardware platforms appear frequently, a modern simulator must be portable, li- 
bRoadRunner has no run-time dependencies beyond standard system libraries and it supports any processor 
LLVM supports. libRoadRunner is written in C++, so it interfaces easily with other C++-bas ed sof t ware, li- 


bRoa dRunner also provides a C language wrapper for cross-language support, and uses SWIG (jBeazlev et al 


1996h to provide a customized native-Python API. The use of SWIG will allow future support for additional 


native language bindings, such as JavaScript, R, or MATLAB, depending on demand. 


Extensibility libRoadRunner’s modular design is easy to maintain and extend. All top-level components, 
such as solvers and integrators, interact via well-defined boundaries (pure virtual interfaces) to reduce inter- 
component dependencies and hide their internal details. A new solver only needs to implement a standard 
interface to function as part of the library, so adding a solver requires no modification of pre-existing code. 


Systems Biology Markup Language as a Declarative Language 

SBML dHucka et al. I . l2003h is a declarative specification format for network models. Because of its history, 
SBML terminology derives from biochemistry and includes common biochemical-reaction abstractions like 
reaction steps, compartments and reaction rate laws, though it can describe any model of form: 


—x(t) = /(x(t),p), 


(1) 


where x is the state vector of the model, and p is a vector of time-independent parameters. 

SBML-specified models can also include events , discontinuous state changes, which trigger under specified 
conditions. libRoadRunner correctly handles SBML-specified events and extends the SBML specification by 
allowing an SBML event to call an arbitrary user-defined function. 

Declarative specification languages, like SBML, define component objects and their interactions, rather 
than defining procedural control flow (*.e., the sequence in which computational operations proceed on 
execution). An SBML specification lists only the network component objects, their interactions and rate 
relations and events which change these interactions and rates, all of which are intrinsic abstractions in 
SBML. Thus, an author writing a model specification in SBML can focus on the underlying biology or 
chemistry of the model rather than on how to implement the model as a simulation. Since SBML does not 
specify the computational operations to implement a model, the control flow, the solvers to use, or how 
to store the model’s elements, an SBML compiler or interpreter must generate them appropriately from 
the SBML specification. Thus, compiling an SBML model specification is more complex than compiling a 
functionally-equivalent model specification in a procedural language. 

SBML model specifications are easier to share than procedural specifications of equivalent models be¬ 
cause they are not implementation dependent; any of the numerous SBML compliant tools can process any 
SBML model speci fi cati on. This porta bility allows model archiving ( e.g ., in exchange repositories such as 
BioModels ( BioModels.net TeamLl2014lf ) and reuse and the relatively simple assembly of multiple SBML- 
specified submodels into larger models. It also simplifies the scientific validation of SBML-specified models 
and ensures that SBML-specified models remain usable, even if the specific software tools that generated 
them fall out of use. 


6 Architecture 

libRoadRunner is a self-contained, easily embedded library with an object-oriented API natively accessible 
in C, C++ and Python (SWIG allows easy extension to other languages). libRoadRunner’s component- 
oriented design specifies a small number of standardized software interfaces (protocols) and how they in¬ 
teract, implemented using standard C++ data types. Component-orientation separates the implementation 
of a component from its interface, so components are easy to add or replace and component swapping 
requires no changes to existing code. E.g., we can add new integrators, steady-state solvers or SBML 
compilers to the libRoadRunner library via the Integrator, SteadyStateSolver and ExecutableModel 
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interfaces respectively. libRoadRunner includes three implementations of the Integrator interface: two 


deterministic integrators (one based on the CVODE integrator from the Sundials suite (jHindmarsh et al . 


2009 ) and the other a standard fourth-order Runge-Kutta method) and a standard Gillespie SSA stochas¬ 
tic inte grator. libR oadR unner implements the SteadyStateSolver interface as a class which uses the 
NLEQ (Nowak and Weim ann. 119911) solver, and we are currently developing additional methods. libRoad- 
Runner implements the ExecutableModel interface as a class which uses our SBML-to-CPU JIT compiler 
(see §0, and we are developing a n SBML-to-GPU J I T com piler as wel l. libRoadRunn e r stat ically links 
to the third-party libraries LLVM ( Lattner and Adve . 2004 ). libSBML ( Bornstein et al. . l2008h . CVODE, 
NLEQ2, LAPACK and Poco. 


7 SBML-to-CPU-Executable Compilation 

LibRoadRunner’s SBML JIT compiler compiles SBML models in the form of s trings to execu t able n ative 
machine code, in memory. Compilation follows the canonical compiler phases ( Alfred V. Aho . fl986h : (1) 
lexical analysis, (2) syntactic analysis, (3) semantic analysis, (4) intermediate code generation, (5) code 
optimization, and (6) native code generation. Standard generic libraries can perform phases 1, 2, 5 and 6. 
However, semantic analysis (phase 3) is specific to the source language. 

In phases 1 and 2, the compiler reads the source text, parses it, and extracts and converts the text’s 
syntactic information into an abstract syntax tree (AST) data structure. Each node in the AST is an essential 
construct such as a n operator, symbo l, literal or function call. Most SBML simulators use components of 
the libSBML ( Bornstein et all 2008 ) library to perform lexical and syntactic analyses of SBML model 
specifications. 

In phases 3 and 4, the compiler reads the AST and assembles it into a sequence of intermediate language 
(IL, a machine-independent assembly language) instructions which form a procedural instantiation of the 
SBML model specification. 

CPUs cannot execute IL programs directly, so phases 5 and 6 optimize the IL (by removing redundant 
operations, optimizing memory layout,...) and convert it into executable machine code. libRoadRunner uses 
components of the LLVM library for phases 5 and 6. An advantage of the LLVM code-generator is that it 
can produce executables for x86, GPUs, ARMs and other processors. 

After the completion of phases 1-6, the JIT compiler returns the executable code in the form of a list of 
callable functions to the calling program. 

During phase 3 (semantic analysis), the compiler must map language symbols to memory address loca¬ 
tions. The compiler of a procedural language such as C, allocates a memory location to each symbol (e.g., 
a variable or function declaration), and resolves that symbol to that location whenever the source code ref¬ 
erences that symbol. Procedural-language compilers map symbols to memory locations using a symbol table 
data structure. SBML has no construct for creating new variables or eliminating variables at run-time, so 
the compiler can compute the exact memory requirements for all symbols and store the symbols in a con¬ 
tiguous memory block. At run-time, during a time-series computation, the libRoadRunner library connects 
a JIT-compiled function to an integrator, which, in turn, calls a function which calculates the rate of change 
of the state vector. Since both the state vector and the rate of change occupy contiguous memory blocks and 
have the same layout as the SBML model variables, the calls pass only two pointers and require no memory 
copying or rearrangement. 

However, compilation of SBML poses challenges. SBML model specifications may define rules which 
state that an expression should replace a specified symbol, or a rate rule which specifies a rate of change 
of the value of a symbol, rather than the symbol value itself. SBML also allows different rules to apply in 
different contexts, such as special rules which only apply when the model is loaded (initial assignment rules). 
Mapping symbol names to memory locations is not one-to-one so a symbol table is insufficient to store the 
mapping. 

Some SBML model simulators allocate storage space for both normal and rule-defined symbols and use 
auxiliary functions to evaluate the rules at run-time as the symbols are read. However, this approach wastes 
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memory storing symbols which resolve to other symbols and complicates execution, as the run-time must 
keep track of rule dependencies. 

Our solution is to extend the symbol table into a symbol forest , a hash table which maps symbol names 
to ASTs describing all the symbols’ rules. The SBML compiler uses the symbol forest much as a procedural- 
language compiler uses a symbol table, to resolve symbol names to memory locations. However, the symbol 
forest must apply any rules which relate symbols to determine the memory location for a given symbol. 
E.g ., if the symbol x has the assignment rule x —> y + 1, whenever the compiler references symbol x , the 
symbol forest will find the rule, generate a sequence of IL instructions which both implement the right hand 
side ( RHS ) of the rule and create a temporary variable to store the result of the rule calculation. The 
symbol forest then stores this sequence of IL instructions and returns the memory location of the instruction 
sequence to the compiler. Later in compilation, the LLVM code generator translates these IL instructions 
into an executable, which calculates and returns the value of the symbol at run-time. The symbol forest 
resolves automatically recursive rules in which the symbols in the RHS of a rule depend on other rules. 

Naively generating IL expansions of the rule definitions inline and creating temporary variables for 
rule evaluation would generate redundant instructions which would slow both compilation and execution. 
libRoadRunner’s scoped symbol cache reduces such redundancy. Many functions in libRoadRunner do not 
modify SBML-defined parameters and variables during function execution, so any rules depending on these 
parameters and variables need evaluation only once during a given call to these functions. Even if the rule 
involves a condition, e..g., x —>• {b if (a > 1) else c), if the function does not change the values of a, b and c, 
the function need evaluate the rule to obtain the value of x only once per call. The SBML compiler therefore 
generates code which evaluates the rule whenever the function is called and stores the result in a temporary 
variable. During a call to the function, the first reference to the symbol evaluates the rule and caches its 
result, and any subsequent references to that symbol during the function call reference the cached value. 
Using a scoped symbol cache reduces both memory usage and execution time, typically by a factor of 10 for 
large models. 

When JIT-compiled functions contain conditional branches which contain rules, the SBML compiler 
generates redundant IL code, which slows compilation (which scales as the size of the IL code) but has no 
speed cost at execution. If the compiler examined all possible branches, determined what rules were present, 
and created temporary variables to contain the results of the rule evaluations, it would reduce the size of 
the resulting IL code, speeding compilation. However, slower execution would offset the faster compilation, 
since the executable would evaluate all rules in all branches, not only those which it needed. We may add a 
compiler directive to allow the user to choose the second option in a future release of libRoadRunner. 


8 Results 


Performance 


Simulation engines which interpret SBML models ( Romer et al. . 1996h . are inherently slower, sometimes 


much slower, than engines which genearte and execute complied code. libRoadrunner uses JIT compilation 
to generate particularly fast simulations. 

Simulation speed depends on the performance of both the state-vector rate calculation and the numeric 
integrator. Since we cannot separate these calculations in most SBML-model packages, we compared an 
SBML model JIT-compiled using libRoadRunner with a hard-coded C-j—(- version of the same model. The 
model implemented 1000 instances of a unimolecular reaction, in which a single substrate goes to a single 
product (S —» P) at a rate of 
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This equation m odel s a reversible H ill equation ( Sauro . 2012 ) that was 
Bowden (IHofmevr and Cornish-BowdenL Il997l ). This test measured the 
between 1000 to 15,000 time steps, and on average, the JIT compiled 
that of a hard-coded C++ specification. All tests were performed on a 
used as the C++ compiler. The timing data implies that we are close to 
than using a completely different approach, such as parallelization, it is 
further. 


developed by Hofmeyr and Cornish- 
total simulation time for performing 
SBML model performance was 96% 
64 bit Linux system, and clang was 
optimal performance. That is, other 
not possible to improve performance 


Python Bindings 

The Python API employs a simple, concise object model, and follows the style and conventions of the widely- 
used SciPy library which makes it easy to learn. The API also provides high performance, low-overhead 
access to the libRoadRunner library. The API only communicates using standard Python data types such 
as lists, dictionaries and NurnPy arrays, which simplifies integration with existing applications. The NumPy 
array type is a data structure which wraps a Python interface around a standard C numeric array. Even 
large NumPy arrays have low overhead, since they return only pointers to internal arrays owned by the 
libRoadRunner library, with no copying of memory. 

libRoadRunner extends the NumPy array to contain row and column name information and to support 
access to rows and columns by name, and extends the NumPy print method to format the extended arrays 
for printing. Since libRoadRunner users frequently need to display the components and interaction names 
in the stoichiometry matrix, the libRoadRunner API makes this display is a single line task: 
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Running a simulation is straight forward and only requires users to load a model and call the simulation 
method. Defaults are in place that include time spans and number of points generated during a simulation. 
By default the simulate method will return time in the first column and all floating model species in the 
remaining columns: 


r = RoadRunner("glycolysis.sbml") 
m = r.simulate(plot=True) 

Here m is a NumPy array, and the optional plot=True argument to the simulate method calls the standard 
plotting library, matplotlib, to display a basic time-series plot of the simulation result. We can customize the 
simulation parameters using optional arguments, e.g., to generate a 100 data-point time series for parameter 
“PI” and concentration “SI” from an SBML-specified model between time t = 0 and time t = 12, we specify: 


r = RoadRunner("glycolysis.sbml") 

m = r.simulate(0, 10, 100, [’time’, ’P’, ’[SI]’]) 


A variety of other built-in symbols are provided to access information on reaction rates, rates of change, 
eigenvalues etc. Like a MATLAB top-level function, the libRoadRunner simulate method provides a con¬ 
sistent front end to all libRoadRunner’s integration engines. Since MATLAB is familiar to many scientists, 
a MATLAB-like architecture should reduce the effort to learn the libRoadRunner API. To simplify genera¬ 
tion of simulation documentation, libRoadRunner meth ods su pport internal pydoc strings which interactive 
Python environments such as IPython or Tellurium ( Sauro et all . 1201 41). make available as pop-up hints. 

In analogy with MATLAB, the libRoadRunner API uses dynamic Python object properties. Loading 
an SBML-specified model via libRoadRunner automatically adds the SBML model’s symbol names to the 
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RoadRunner object, allowing dynamic introspection and modification of the object. If a model contains pa¬ 
rameters and species ’ x ’, ’y ’, ’ SI 1 , ’ S2 ’, the RoadRunner object will include these names as properties, 
which a user can read or set. E.g., 


# load a model that has ids ’x’, ’y’ 

r = RoadRunner(’some.model.xml’) 
r.x = 1.5 # set the ’x’ parameter 

r.y = 2.0 # set the ’y’ parameter 

print(r.Sl) # print the ’SI’ species 


and ’SI ’ 

to 1.5 
to 2. 0 

concentration 


Support for Analysis 


Since libRoadRunner was inspired and developed from the original C# roadrunner, libRoadRunner inherits 
many of the analysis functions that roadRunner has. This includes methods to calculate both scaled and 
unsealed control coefficients, elasticities, rates to changes in all parameters, including conserved qu antities , 
eigen values and eigenvectors and a variety of stoichiometric results such as the Link and K matrices (Reder, 
19881) . libRoadRunner can also compute the frequency response in order to make Bode plots. 


Identification of Conserved Quantities 

Many computations of biochemical networks requires identification of conserved quantities (moieties in bio- 
chemic al usage) a nd e limination of linearly-dependent species to avoid inversion of singular Jacobian ma¬ 


trices ( Vallabhaiosvula et al . 20061) . libRoadRunner implements a libSBML plug-in which performs this 
reduction on a SBML Document object, first identifying conserved quantities and dependent species, then 
adding the conserved quantities to the document as set of global parameters and replacing the dependent 
species with assignment rules. The user can modify these conserved quantities, which behave as parameters, 
to investigate their effect on the dynamics of the model. 


Hybrid Models 

libRoadRunner includ es thr ee independent simulation approaches, deterministic, stochastic and Flux Balance 
Analysis ( Orth et "oil l20icl) . As a result lib RoadR unne r is ideally suited for simulating hybrid models such 


as the whole-cell hybrid model by Karr (jKarr et al. 
through python connecting code. 


2012) where the different methods can be combined 


9 Use Cases 

libRoadRunner’s ease of use, ability to handle very complex SBML models and very fast model execution 
speed have led to its rapid adoption in a variety of applications. 


Interactive Application 


A cross-platform integrate d Python environment called Tellurium (jSauro et al. 120141) has been developed 


based on the Spyder IDE ( Ravbaut and Cordoba . 20151) . which combines libRoadRunner, libSBML, Anti¬ 
mony ( Smith et «Ll . l2009lh libSEDML and a variety of other packages. Tellurium provides a comprehensive 
development and analysis environment for SBML-specified models. Tellurium uses Antimony as a model 
specification language and libRoadRunner provides simulation execution and analysis functions. While li- 
bRoadRunner’s execution speed is not a priority in this application, the concise syntax and intuitiveness of 
its Python API are essential to Tellurium’s convenient interactive creation, simulation and analysis of DMN 
models. A detailed description of Tellurium will be given in a separate publication. 
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Integrating SBML-model specifications into Multicell Models simulated in Com- 
puCell3D 

CompuCell3D ( CC3D ), a simulation environment for multiscale, multicell virtual-tissue model development 
and simulation, was the first tool to adopt libRoadRunner as a core engine. CC3D proper defines an object 
class of cell agent and behavior methods to allow cell agents to grow, divide, die, secrete/absorb chemicals, 
move, etc.... The libRoadRunner/CC3D integration allows the state of an SBML-specified model inside a 
cell object to control the CC3D parameters describing the cell agent’s behaviors and vice versa. CC3D also 
uses libRoadRunner to time-step the state of a biochemical models and handles their interactions with other 
model objects via Python. 

E.g ., in a model of changes in cell-cell adhesion leading to invasive tumor phenotypes, the CC3D cell 
objects have a CC3D parameter adhesion-molecule density , which controls the CC3D behavior cell-cell adhe¬ 
sion. An SBML-specified model rel ates t he level o f t he tr ansmembrane adhesion receptor E-cadherin in each 
cell to the cells’ level of /3-catenin ( Andasari et all 1 201 2l) . The CC3D-model specification uses the libRoad¬ 
Runner Python API to connect the CC3D adhesion-molecule density to the SBML-model transmembrane 
E-Cadherin level. At run-time, libRoadRunner engine time evolves the model inside the cells, while the 
CC3D engine handles the evolution of the cell objects. 

Another example of the use of SBML models in virtual-tissue modeling is simulation of Delta-Notch 
patterning during embryonic development. Delta and Notch are heterophilic transmembrane receptors whose 
signaling is mutually inhibitory within a cell. However, the level of signaling depends on both the amount 
of Delta on the membrane of a cell and the amount of Notch on the surfaces of neighboring cells and 
vice versa. Thus, to understand the dynamics of the signaling network, we need to know not only the 
model within the cell, but the pattern of its contacts with neighboring cells and their levels of Delta and 
Notch. To model this situation, we create CC3D cell objects and arrange them in an epithelium (a quasi-2D 
sheet). Each cell contains an SBML-specified model that describes how the cell’s levels of membrane-bound 
and cytosolic De lta and No tch change, given a particular input level of transmembrane Delta and Notch 


signaling ([Swat et all, 2012). A Python layer uses the libRoadRunner API to combine the amount of Delta 


on the membrane of each cell, the amount of Notch on the membrane of each adjacent cell (adjacency is a 
CC3D model parameter) and the (CC3D model) contact area between each adjoining pair of cells to calculate 
the strength of Delta and Notch signaling each cell experiences. libRoadRunner then updates the internal 
state of the Delta-Notch signaling and regulatory networks using these signaling strengths as boundary 
conditions, while CC3D updates the cellshapes and positions and identifies cell rearrangements changing 
adjacencies and contact areas between cells. Together, the interactions produce the typical checkerboard 
pattern of embryonic Delta-Notch signaling. 


Multi-scale Modeling of Liver Metabolism using libRoadRunner for Numerical 
Integration 

Liver function arises as a complex interaction among morphology, perfusion, and metabolism on multiple 
scales. E.g., the Virtual Liver Network has developed an organ-level model of human galactose clearance 
which includes single cell metabolism of hepatocytes, a representation of ultra-structure and micro-circulation 
in hepatic tissue, and a description of the structure of the entire organ to investigate how galactose cle aran ce 
on organ level results from the interplay of blood flow, tissue structure and cellular metabolism (Konig, 
20151b 

The liver model includes an SBML-specified model of the smallest functional unit of the liver, a sinusoid 
consisting of a perfused capillary surrounded by hepatocytes. The sinusoid model contains a biochemical 
network describing galactose metabolism in individual hepatocytes, coupled via SBML-specified discretized 
transport equations for convection and diffusion, resulting in a model with several thousand components 
and interactions. The model also uses SBML events to describe the time-varying supply of galactose to 
the liver. To account for heterogeneity in blood flow and tissue architecture required simulation of more 
than 10 5 (around 100,000 - 200,000 simulations necessary in total) replicas of the model under varying 
tissue and flow parameters, which became feasible only due to libRoadRunner’s fast time-series generation 












and support for variable step sizes, which dramatically reduced output file size. libRoadRunner’s Python 
API supported rewrite-free integration of the SBML models into a complex pre-existing modeling workflow, 
which included data management using Django, model annotations using Python bindings to libSBML, 
model prototyping using Python bindings to antimony and visualization of results using the Python REST 
interface to Cytoscape (CySBML, CyFluxViz). 


Modeling of Synaptic, Neuronal and Neuron Network Dynamics in MEMORY 
platform using libRoadRunner 


The M EMORY platform (Multiscale intEgrated Model Of the neRvous sYstern, formerly EONS (I Bouteiller et al 


120081) 1 simulates the function and dynamics of elements ranging from single channels or receptors (elementary 
models), to synapses, which include many elementary models, to neurons, which themselves may include a 
large number of synapses. E.g., an SBML-specified neuron model may include many SBML-specified synapse 
models, each of which includes many SBML-specified neurotransmitter release and diffusion, AMPA recep¬ 
tor and NMDA-receptor models (both ionotropic receptors for the glutamate neurotransmitter). MEMORY 
depends on the flexibility and ease of use of libRoadRunner to enable the assembly of such complex hier¬ 
archical inodels^Neuronal models may be large, e.g., representing 10 ionotropic synapses in a CA1 neuron 
model ( Izhikevichl . 120031) requires 73 events, 290 reactions, 414 rules and 1459 parameters. The speed of 
libRoadRunner time-series generation enables MEMORY to solve complex neuronal models quickly. 

Additionally, to ensure that a neuronal model quantitatively predicts biological functions like mem¬ 
brane potentials or intracellular molecular concentrations, MEMORY can optimize the model’s parameters 
by fitting between multiple simulation and experimental time-series for characteristics including changes 
in receptor conductance, desensitization properties and spiking patte rns. MEMORY uses evolutionary 


multi-objective optimization (from the EMOO framework (Bahl et al J, [2012)), which requires large num¬ 
bers of simulation replicas. E.g., elementary-model optimization of an NMDA-receptor model with re¬ 
spect to eight distinct experimental results for dynamical changes in receptor channel conductance following 
paired-pulse stimulation, required 15,000 generations with 400 individuals per generation, i.e., 6 million 
simulation replicas (corresponding to 13,000 hours of simulated time). libRoadRunner took 66 hours to 
run the ent ire optimization on a 400-node computer cluster, orders of magnitude faster than other SBML 
simulators ( Bouteiller et <aZ.I . l2015h . 


10 Conclusions 

libRoadRunner’s speed and ease of integration allow researchers to solve very large models, to include models 
in multiscale systems and run large ensembles of smaller models. The use of Python makes simulations 
easy to carry out for new users, while C++ and C APIs are attractive to developers looking to integrate 
libRoadRunner’s capabilities into their existing simulation frameworks. libRoadRunner runs on both x86 
and ARM architectures. We have successfully deployed libRoadrunner to Windows, Mac OS X, Linux, 
Raspberry Pi, NVIDIA Jetson TK1 and Parallella boards. 


11 Future Work 

Our roadmap currently includes the following proposed extensions to libRoadRunner. 


Support Automatic Differentiation of Jacobian calculations libRoadRunner’s stability and time 
evolution modules require numerical differentiation of the evolution equations to calculate a Jacobian ma¬ 
trix. For stiff equations, numerical differentiation requires the use of slow implicit solvers, which can introduce 
round-off errors. While we could sometimes calculate the Jacobia n analytically from the AST, symbolic dif¬ 
ferentiation can produce equations with unpractically many terms ( San 10 . 1993 ). We are therefore evaluating 
an intermediate solution, known as automatic differentiation (AD) ( Griewank et al. . 19891) . which computes 
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analytic derivatives of programs rather than mathematical expressions. Simple analytical derivatives ex¬ 
ist for the elementary mathematical functions (such as sin, cos, exp, etc.) available in most programming 
languages. Since most functions a program computes are compositions of these intrinsic functions, their 
derivatives are chain-combinations of the elementary derivatives. A JIT compilation of the resulting pro¬ 
gram for the Jacobian is usually much more compact than the equivalent analytical expansion. 


Improve Steady-State Solvers libRoadRunner uses the FORTRAN NLEQ2 nonlinear steady-state 
solver, which is not thread safe (multiple execution threads may call a routine concurrently). The re¬ 
mainder of the libRoadRunner library is thread-safe, so we had to place exclusive access locks (mutexes) on 
the NLEQ solver which restricts its use to one thread at a time. To elimin ate thi s restricti on, we plan to 
add several thread-safe ste ady state solvers, including Sund ials’ KINSOL ( Hindmarsh et all l2005lh and a 


damped Newton algorithm ( Hoops et al . 20061 Saurol . 1993 b 


Add Additional Integrators libRoadRunner includes the explicit and implicit integrators from the 
Sundials suite, a Ru nge-Kutta integrator and a Gillespie integrator. We plan to add integrators based on 
LSODA, ARKode ( Reynolds et al\ . 2014 ) for the Sundials suite, and the banded integrator from CVODE 
(plus a tool to detect Jacobians with banded structure). The CVODE integrator is much faster and uses 
much less memory than other integrators when solving models with banded Jacobians. 


JDesigner A new cros splafor m version of JDesigner is under development and will employ libRoadRunner 
as its simulation engine (i Saurol 20151) . 
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